/*   AUTHOR: Zachary D. Taylor, Stephen Durfey, Matthew Wehinger

     DATE: 11/11/2010

     PURPOSE: This program uses openmp to perform matrix multiplication on arrays of 
	      integers, doubles, and floating point variables. The execution times are 
	      compared to each other and to a similar program written using CUDA.

     command line arguments: a.out thread_number m n k
*/

#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <omp.h>
#include <iomanip>

using std::cin;
using std::cout;
using std::endl;
using std::right;
using std::setw;

//  Function prototypes
template <class T>
void initializeMatrix(T **, int, int);
template <class T>
void populateMatrix(T **, int, int);
template <class T>
void matrixMultiply(T **, T **, T **, int, int, int);
template <class T>
void deleteMatrix(T **, int);
double calculateSpeedup(double, double);
template <class T>
void printMatrix(T **, int, int);

// Begin main function
int main(int argc, char** argv) {
	//  Matrix dimensions
	int m = 0, n = 0, k = 0;

	m = atoi(argv[2]);
	n = atoi(argv[3]);
	k = atoi(argv[4]);

	//  Time related variables
	double startTime, endTime;
	double executionTimeI, executionTimeD, executionTimeFP;

	//  Dynamic arrays
	int ** matrixA, ** matrixB, ** matrixC;

	//  Allocate space for arrays
	matrixA = new int*[m];
	matrixB = new int*[n];
	matrixC = new int*[m];

	if(!matrixA || !matrixB || !matrixC) {
		cout << "Memory allocation failed.  Exiting..." << endl;
		return(1);
	}

	//  Initialize matrices to zero
	cout << "Initializing matrices..." << endl;
	initializeMatrix(matrixA, m, n);
	initializeMatrix(matrixB, n, k);
	initializeMatrix(matrixC, m, k);

	//  Populate matrices with random values
	cout << "Populating matrices..." << endl;
	populateMatrix(matrixA, m, n);
	populateMatrix(matrixB, n, k);

	//  Multiply matrices in parallel
	//  Integer
	//  Measure execution time
	cout << "Multiplying integer array..." << endl << endl;
	startTime = omp_get_wtime();
	matrixMultiply(matrixA, matrixB, matrixC, m, n, k);
	endTime = omp_get_wtime();
	executionTimeI = endTime - startTime;
	
	//  Free memory
	deleteMatrix(matrixA, m);
	deleteMatrix(matrixB, n);
	deleteMatrix(matrixC, m);

	//  Dynamic arrays
	double ** matrixD, ** matrixE, ** matrixF;

	//  Allocate space for arrays
	matrixD = new double*[m];
	matrixE = new double*[n];
	matrixF = new double*[m];

	if(!matrixD || !matrixE || !matrixF) {
		cout << "Memory allocation failed.  Exiting..." << endl;
		return(1);
	}

	//  Initialize matrices to zero
	cout << "Initializing matrices..." << endl;
	initializeMatrix(matrixD, m, n);
	initializeMatrix(matrixE, n, k);
	initializeMatrix(matrixF, m, k);

	//  Populate matrices with random values
	cout << "Populating matrices..." << endl;
	populateMatrix(matrixD, m, n);
	populateMatrix(matrixE, n, k);

	//  Multiply matrices in parallel
	//  Double
	//  Measure execution time
	cout << "Multiplying double array..." << endl << endl;
	startTime = omp_get_wtime();
	matrixMultiply(matrixD, matrixE, matrixF, m, n, k);
	endTime = omp_get_wtime();
	executionTimeD = endTime - startTime;
	
	//  Free memory
	deleteMatrix(matrixD, m);
	deleteMatrix(matrixE, n);
	deleteMatrix(matrixF, m);

	//  Dynamic arrays
	float ** matrixG, ** matrixH, ** matrixI;

	//  Allocate space for arrays
	matrixG = new float*[m];
	matrixH = new float*[n];
	matrixI = new float*[m];

	if(!matrixG || !matrixH || !matrixI) {
		cout << "Memory allocation failed.  Exiting..." << endl;
		return(1);
	}

	//  Initialize matrices to zero
	cout << "Initializing matrices..." << endl;
	initializeMatrix(matrixG, m, n);
	initializeMatrix(matrixH, n, k);
	initializeMatrix(matrixI, m, k);

	//  Populate matrices with random values
	cout << "Populating matrices..." << endl;
	populateMatrix(matrixG, m, n);
	populateMatrix(matrixH, n, k);

	//  Multiply matrices in parallel
	//  Floating Point
	//  Measure execution time
	cout << "Multiplying floating point array..." << endl << endl;
	startTime = omp_get_wtime();
	matrixMultiply(matrixG, matrixH, matrixI, m, n, k);
	endTime = omp_get_wtime();
	executionTimeFP = endTime - startTime;
	
	//  Free memory
	deleteMatrix(matrixG, m);
	deleteMatrix(matrixH, n);
	deleteMatrix(matrixI, m);

	cout << "Done." << endl << endl;

	//  Print results
	cout << "Execution Time" << endl;
	cout << "--------------" << endl;
	cout << "Integer Multiplication:                  " << executionTimeI << " seconds" << endl;
	cout << "Double Multiplication:                   " << executionTimeD << " seconds" << endl;
	cout << "Floating Point Multiplication:           " << executionTimeFP << " seconds" << endl << endl;

	cout << "Speedup" << endl;
	cout << "--------------" << endl;
	cout << "Double:                   " << calculateSpeedup(executionTimeI, executionTimeD) << endl;
	cout << "Floating Point:           " << calculateSpeedup(executionTimeI, executionTimeFP) << endl << endl;

	system("PAUSE");

	return 0;
}

//  Initializes matrix of size numRows x numCols
//  Values set to zero
template <class T>
void initializeMatrix(T **matrix, int numRows, int numCols) {
	for(int i = 0; i < numRows; i++)
		matrix[i] = new T[numCols];

	for(int i = 0; i < numRows; i++)
		for(int j = 0; j < numCols; j++)
			matrix[i][j] = (T) 0;
}

//  Populates matrix of size numRows x numCols
//  Values are floating point values between -10 and 10
template <class T>
void populateMatrix(T **matrix, int numRows, int numCols) {
	for(int i = 0; i < numRows; i++)
		for(int j = 0; j < numCols; j++)
			matrix[i][j] = (T) -10 + rand() * 20 / RAND_MAX;
}

//  Multiplies two matrices of size m x n and n x k
//  Stores result in matrix of size m x k
//  Pragma tells compiler that program should execute in parallel
//  Uses static scheduling/partitioning
template <class T>
void matrixMultiply(T **matrixA, T **matrixB, T **matrixC, int m, int n, int k) {
	int chunkSize = 0, chunkRemainder = 0;

	#pragma omp parallel
	{
		chunkSize = m / omp_get_num_threads();
		chunkRemainder = m % omp_get_num_threads();

		int id = omp_get_thread_num();
		int start = id * chunkSize, end = start + chunkSize;

		for(int p = start; p < end; p++)
			for(int q = 0; q < n; q++)
				for(int r = 0; r < k; r++)
					matrixC[p][r] += matrixA[p][q] * matrixB[q][r];
	}

	if(chunkRemainder > 0) {
		for(int p = m - 1; p >= m - chunkRemainder; p--)
			for(int q = 0; q < n; q++)
				for(int r = 0; r < k; r++)
					matrixC[p][r] += matrixA[p][q] * matrixB[q][r];
	}
}

//  Frees memory allocated to matrix with numRows rows
template <class T>
void deleteMatrix(T **matrix, int numRows)	{
	for(int i = 0; i < numRows; i++) {
		delete[] matrix[i];
		matrix[i] = NULL;
	}

	delete[] matrix;
	matrix = NULL;
}

//  Calculates speedup for the given execution times
double calculateSpeedup(double executionTime1, double executionTime2) {

	if(executionTime2 == 0)
		return 1;
	else
		return executionTime1/executionTime2;
}

template <class T>
void printMatrix(T **matrix, int numRows, int numCols) {
	for(int i = 0; i < numRows; i++) {
		for(int j = 0; j < numCols; j++) {
			cout << setw(5) << right;
			cout << matrix[i][j];
		}
		cout << endl;
	}
	cout << endl;
}
